home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 06 - 1990 / 06.11 Nov 90 / NeuralNet Estimator⁄EDIT / Gradient methods.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-03  |  14.2 KB  |  436 lines  |  [TEXT/KAHL]
  1. #include "Neural Network.h"
  2. #include <math.h>
  3.  
  4. extern FILE * Jac;
  5.  
  6. extern NeuralNet *    theNet;
  7. extern DTypeVector     yData;
  8. extern DTypeMatrix    XData;
  9. extern DTypeVector     Alpha[];
  10. extern DTypeMatrix  Jac_T;
  11. extern DTypeVector    Pi;
  12. extern DTypeVector    Diag;
  13. extern DTypeVector     Resid;
  14. extern DTypeVector     dSquash[];
  15. extern DTypeVector    gradSS;
  16. extern DTypeMatrix  Phi, T2;
  17.  
  18. static unsigned int    totparms;    /* total parms in model, set by AllotGradientWorkSpace() */
  19.  
  20. /*-------------------------
  21.     Find parameter values that minimize sum of squares, uses line search method.
  22.     Expects:
  23.         1. tolerance values gradtol, steptol; which are stored in NeuralNet structure.
  24.         2. validly defined and alloted NeuralNet structure.
  25.         3. input data in matrix XData and output data values in vector yData.
  26.     Returns termcode = FAIL if current parm value is not an approximate critical point 
  27.                      = METGRADTOL if norm of scaled gradient less than gradtol
  28.                      = METSTEPTOL if scaled step is less than steptol
  29.                      = LINEFAIL if linesearch failed to find next parm distinct from current value
  30.                      = EXCEDITNLIM if iteratation limit exceeded
  31.                      = SINGJAC if Jacobian is singular
  32. */
  33. do_GaussNewton()
  34. {
  35.     int itncount = 0;    /* start iteration count at zero */
  36.     int termcode = FAIL;/* termcode must be changed to nonzero in order to terminate */
  37.     int retcode;
  38.     int sing;            /* flag for singularity of Jacobian */
  39.     double SS = 0;
  40.     
  41.     while(termcode == FAIL)
  42.     {    itncount += 1;                            /* increment iteration count */
  43.         SS = Compute_SS_Jac();                /* this also computes the Jacobian, stored as transpose */
  44.         Matrix_by_Vec(&Jac_T, &Resid, &gradSS); 
  45.                 /* Compute the gradient of the SS function = Jac_T*Resid. 
  46.                     should be before next step calculation since that computation
  47.                     affects the values in Resid    */
  48.         sing = OLSbyQRmethod(&Jac_T,&Pi,&Diag,&Resid);    
  49.                 /* Compute next step, step increment stored in Diag sing TRUE => singular
  50.                     Jacobian, so stop    */
  51.                     
  52.         if(sing)
  53.             termcode = SINGJAC;
  54.         else
  55.         {    SaveParms(&Pi);        /* save weights in case linesearch fails,the Pi 
  56.                                         vector is not being used at this point, so is 
  57.                                         reused to save space */
  58.             if(theNet->usemaxstep == TRUE)
  59.                 ConstrainStep();            /* constrain step size to maxstep */
  60.             retcode = LineSearch(&SS);    /* also sets the SS value */
  61.             termcode = StopYet(SS,retcode,itncount);
  62.         }
  63.     }/* end of while(termcode == FAIL) */
  64.     return(termcode);
  65. }
  66.  
  67. /*-------------------------
  68.     Find parameter values that minimize sum of squares, uses quasi Gauss-Newton method 
  69.     combined with line search.    Adds a nonsingular diagonal matrix to Jacobian to overcome
  70.     singularity problems.  Takes more memory than Gauss-Newton.  Similar to Levenberg-
  71.     Marquardt but trust region is a constant since only trying to fix singularity of
  72.     Jacobian problem.
  73.     Expects:
  74.         1. tolerance values gradtol, steptol; which are stored in NeuralNet structure.
  75.         2. validly defined and alloted NeuralNet structure.
  76.         3. input data in matrix XData and output data values in vector yData.
  77.         4. total number of parameters for model in "totparms".
  78.     Returns termcode = FAIL if current parm value is not an approximate critical point 
  79.                      = METGRADTOL if norm of scaled gradient less than gradtol
  80.                      = METSTEPTOL if scaled step is less than steptol
  81.                      = LINEFAIL if linesearch failed to find next parm distinct from current value
  82.                      = EXCEDITNLIM if iteratation limit exceeded
  83.                      = NETSAT if network is possibly oversaturated
  84. */
  85. do_quasiGaussNewton()
  86. {    
  87.     int i;
  88.     int itncount = 0;    /* start iteration count at zero */
  89.     int termcode = FAIL;    /* termcode must be changed to nonzero in order to terminate */
  90.     int retcode;
  91.     int sing;            /* flag for singularity of Jacobian */
  92.     double SS = 0;
  93.     
  94.     while(termcode == FAIL)
  95.     {    itncount += 1;                /* increment iteration count */
  96.         ClearDTypeMatrix(&Jac_T);
  97.         SS = Compute_SS_Jac();    /* this also computes the Jacobian, stored as transpose */
  98.         AppendResidZeros();            /* append "totparms" zeros to Resid vector */
  99.         AppendMuIdentity();        /* append Mu times Identity matrix to end of Jacobian */
  100.         ComputegradSSforqGN(); 
  101.                 /* Compute the gradient of the SS function = Jac_T*Resid. Should be before
  102.                     calculation of next step since that computation affects the values in 
  103.                     Resid vector.
  104.                     */
  105.         sing = OLSbyQRmethod(&Jac_T,&Pi,&Diag,&Resid);    
  106.                 /* Compute next step, step increment stored in Diag sing TRUE => singular 
  107.                     Jacobian, so stop */
  108.         if(sing)
  109.             termcode = NETSAT;
  110.         else
  111.         {    SaveParms(&Pi);        /* save weights in case linesearch fails,the Pi 
  112.                                         vector is not being used at this point, so is 
  113.                                         reused to save space */
  114.             if(theNet->usemaxstep == TRUE)
  115.                 ConstrainStep();            /* constrain step size to maxstep */
  116.             retcode = LineSearch(&SS);    /* also sets the SS value */
  117.             termcode = StopYet(SS,retcode,itncount);
  118.         }
  119.     } /* end of while(termcode == FAIL) */
  120.     return(termcode);
  121. }
  122.  
  123. /*-------------------------
  124.     Append zeros to the Resid vector.
  125. */
  126. static
  127. AppendResidZeros()
  128. {
  129.     int i;
  130.     DataType * v;    /* pointer to elements of residual vector */
  131.     
  132.     v = *Resid.cells + XData.rows;    /* XData.rows is the number of observations */
  133.     
  134.     for(i=0; i<totparms; i++, v++)
  135.         *v = 0.0;
  136. }
  137.  
  138. /*-------------------------
  139.     Special function to compute gradient for quasi Gauss-Newton method.
  140. */
  141. static
  142. ComputegradSSforqGN()
  143. {
  144.     register int    i;
  145.     register int    j;
  146.     register int    R;            /* number of rows in matrix    */
  147.     register int    C;            /* number of columns in matrix */
  148.     register DataType temp;        /* ? this uses a floating pt register on 68881 */
  149.     register DataType * y;        /* following declarations use up the 3 address registers */
  150.     register DataType * row;
  151.     register DataType * vec;
  152.  
  153.     R = Jac_T.rows;
  154.     C = Jac_T.cols;
  155.     vec = *Resid.cells;
  156.     y = *gradSS.cells;
  157.     for(i=0; i<R; y++, i++)
  158.     {    row = *Jac_T.cells + i*C;
  159.         temp = (*row)*(*vec);
  160.         for(j=1; j<XData.rows; j++)         /* use XData.rows here instead of Jac_T.cols */
  161.             temp += (*(row+j))*(*(vec+j));
  162.         *y = temp;
  163.     }
  164. }
  165.  
  166. /*-------------------------
  167.     Append a matrix Mu times the Identity matrix to the Jacobian.
  168. */
  169. static
  170. AppendMuIdentity()
  171. {
  172.     register int i,j;
  173.     register DataType * jcell;    /* pointer to cells of Jacobian matrix, stored as transpose */
  174.     register DataType MuI = .001;
  175.         
  176.     /* In the case of quasi Gauss-Newton method the Jacobian has appended to it
  177.         a (totparms x totparms) Identity matrix.  Thus, 
  178.         *Jac_T.cells + XData.rows + i*Jac_T.cols points to the first element of the 
  179.         ith row of this identity matrix, which is the row i, column (XData.row +1) of 
  180.         the Jac_T matrix.
  181.     */
  182.     for(i=0; i<totparms; i++)
  183.     {    jcell = *Jac_T.cells + XData.rows + i*Jac_T.cols;
  184.         for(j=0; j<i; j++, jcell++)
  185.             *jcell = 0.0;
  186.         *jcell = MuI;    
  187.         j++ , jcell++;
  188.         for(; j<totparms; j++, jcell++)
  189.             *jcell = 0.0;
  190.     }
  191. }
  192.  
  193. /*----------------------
  194.     Constrain the step of a Gauss-Newton or quasi Gauss-Newton method.
  195. */
  196. static
  197. ConstrainStep()
  198. {
  199.     register int i;
  200.     register int N = Diag.rows;
  201.     register double snorm = 0.0;
  202.     register double pnorm = 0.0;
  203.     register DataType * cell;
  204.     register double K;
  205.     
  206.     cell = *Diag.cells;
  207.     for(i=0; i<N; i++, cell++)
  208.         snorm += pow((double)*cell,2.0);
  209.     cell = *Pi.cells;
  210.     for(i=0; i<N; i++, cell++)
  211.         pnorm += pow((double)*cell,2.0);
  212.     K = (theNet->maxrelstep*pnorm/snorm) + (theNet->maxstep/sqrt(snorm));
  213.     cell = *Diag.cells;
  214.     for(i=0; i<N; i++, cell++)
  215.         *cell = (*cell)*K;
  216. }
  217.  
  218.  
  219. /*-----------------------
  220.     Alg A6.3.1 of Dennis and Shanabel.
  221.     Expects full Gauss-Newton step in vector Diag.
  222.     Must set the SS value.
  223.     Returns retcode = LINEOK if satisfactory new parm value found
  224.                     = LINENOTOK if can't find a step small enough to reduce SS 
  225. */
  226. static
  227. LineSearch(ss_)
  228.     double * ss_;            /* pointer to value of Sum of Squares */
  229. {
  230.     double ss;                /* value of Sum of Squares */
  231.     double minlambda;        /* min lambda step allowed, smaller than this and would meet stop criteria anyway */
  232.     double initslope;        /* initial slope of SS function */
  233.     double lambda = 1.0;    /* step size, start off with full step */
  234.     int retcode = 2;        /* return code for search */
  235.     double a = .0001;        /* constant used in step acceptance test, called alpha in source,
  236.                                 this value recommended by Dennis and Shanabel */
  237.  
  238.     initslope = -Vec_by_Vec(&gradSS,&Diag);    
  239.                             /* Compute the initial slope of SS function,
  240.                             the Diag function contains the negative of the next step, 
  241.                             so to get the required negative initial slope would need to
  242.                             multiply elements of Diag by -1, instead put a negative sign
  243.                             in front of the calculation 
  244.                             */
  245.     minlambda = Compute_minlambda();        /* Compute the minimum lambda */
  246.     SaveParms(&Pi);            /* the Pi vector is not being used at this point, so is 
  247.                                 reused to save space */
  248.     while(retcode > 1) 
  249.     {    
  250.         UpdateParms(lambda);
  251.         ss = Compute_SS();
  252.         if(ss < *ss_ + a*lambda*initslope)    /* true => satisfactory new parm value found */
  253.         {    retcode = LINEOK;
  254.             *ss_ = ss;            /* reset the Sum of Squares value before exiting */
  255.         }
  256.         else if(lambda < minlambda)
  257.         {    retcode = LINENOTOK;
  258.             RestoreParms(&Pi);    /* restore previous iteration's weight values for network
  259.                                     don't need to reset SS since minimization can't continue */
  260.         }
  261.         else 
  262.             lambda = .1*lambda;
  263.     }
  264.     return(retcode);
  265. }
  266.  
  267. /*-----------------------
  268.     Update the weight parameters by lambda times the step given in vector Diag.
  269.     Expects old weight values to be in the vector Pi.
  270. */
  271. static
  272. UpdateParms(lambda)
  273.     DataType lambda;
  274. {
  275.     int j,k,N;
  276.     DataType * p;    /* pointer to step value */
  277.     DataType * w_;    /* pointer to previous weight values */
  278.     DataType * w;    /* pointer to weight matrix */
  279.     
  280.     p = *Diag.cells;
  281.     w_ = *Pi.cells;
  282.     for(j=0; j<theNet->OutLayer; j++)
  283.     {    N = (theNet->W[j].rows)*(theNet->W[j].cols);
  284.         w = *theNet->W[j].cells;
  285.         for(k=0; k<N; k++, w++, w_++, p++)
  286.             *w = *w_ - lambda*(*p);
  287.     }
  288. }    
  289.     
  290. /*-------------------------
  291.     Determine if should stop searching for minimum.
  292.     Expects step in vector Diag.
  293.     Returns termcode = FAIL if current parm value is not an approximate critical point 
  294.                      = METGRADTOL if norm of scaled gradient less than gradtol
  295.                      = METSTEPTOL if scaled step is less than steptol
  296.                      = LINEFAIL if linesearch failed to find next parm distinct from current value
  297.                      = EXCEDITNLIM if iteratation limit exceeded
  298. */
  299. static
  300. StopYet(ss,retcode,itncount)
  301.     double ss;        /* the Sum of Squares value should be set by LineSearch() */
  302.     int retcode;    /* return code from LineSearch() */
  303.     int itncount;    
  304. {
  305.     int j,k,N;
  306.     DataType * v;            /* pointer to cells in vectors gradSS  or Diag*/
  307.     DataType * w;            /* pointer to cell of weight matrix */
  308.     double rel = 0.0;        /* hold value of relative gradient */
  309.     double max = 0.0;        /* hold maximum relative gradient value */
  310.     int        termcode = FAIL;
  311.     
  312.     if (retcode == LINENOTOK) 
  313.         termcode = termcode | LINEFAIL;
  314.     
  315.     /*---- First check for zero gradient ----*/
  316.     v = *gradSS.cells;
  317.     for(j=0; j<theNet->OutLayer; j++)
  318.     {    N = (theNet->W[j].rows)*(theNet->W[j].cols);
  319.         w = *theNet->W[j].cells;
  320.         for(k=0; k<N; k++, v++, w++)
  321.         {    rel = fabs((double)((*v))*(fabs((double)(*w))/ss));    /* don't need abs(ss) since ss is always positive */
  322.             max = (max < rel) ? rel : max;
  323.         }
  324.     }
  325.     if (max < theNet->gradtol) termcode =(termcode | METGRADTOL);
  326.     
  327.     /*---- Second check for zero step ----*/
  328.     v = *Pi.cells;
  329.     max = 0.0;
  330.     for(j=0; j<theNet->OutLayer; j++)
  331.     {    N = (theNet->W[j].rows)*(theNet->W[j].cols);
  332.         w = *theNet->W[j].cells;
  333.         for(k=0; k<N; k++, v++, w++)
  334.         {    rel = fabs(((double)(*v-*w))/((double)(*w)));    /* don't need abs(ss) since ss is always positive */
  335.             max = (max < rel) ? rel : max;
  336.         }
  337.     }
  338.     if (max < theNet->steptol) termcode =(termcode | METSTEPTOL);
  339.     
  340.     if (itncount > theNet->itnlimit) termcode = (termcode | EXCEDITNLIM);
  341.  
  342.     return(termcode);
  343. }
  344.  
  345.     
  346. /*-----------------------
  347.     Compute the minimum lambda allowed for line search.
  348.     Any lambda value lower than this value would meet the stop criteria for minimum
  349.     step anyway.
  350. */
  351. static double
  352. Compute_minlambda()
  353. {
  354.     int j,k,N;
  355.     DataType * p;            /* pointer to cells in Diag, the proposed step size */
  356.     DataType * w;            /* pointer to cell of weight matrix */
  357.     double rellength = 0.0;    /* hold kth value of relative gradient */
  358.     double maxrel = 0.0;    /* hold maximum relative gradient value */
  359.  
  360.     p = *Diag.cells;
  361.     for(j=0; j<theNet->OutLayer; j++)
  362.     {    N = (theNet->W[j].rows)*(theNet->W[j].cols);
  363.         w = *theNet->W[j].cells;
  364.         for(k=0; k<N; k++, p++, w++)
  365.         {    rellength = fabs((double)((*p)/(*w)));    
  366.             if (maxrel < rellength) maxrel = rellength;
  367.         }
  368.     }
  369.     rellength = theNet->steptol/maxrel;    /* just using rellength to calculate return value */
  370.     return(rellength);
  371. }
  372.  
  373. /*----------------------
  374.     Allot memory for data structures used by the Jacobian matrix and other structures used
  375.     in minimization.
  376.     Physical storage of Jacobian is as the transpose.
  377.     Requires # observations from the data structure.
  378.     Since is always run before method, also sets the totparms variable.
  379. */
  380. AllotGradientWorkSpace()
  381. {
  382.     int i;
  383.     int mxprms = 1;        /* keep track of layer with largest number of parameters */
  384.     
  385.     totparms = 0;        
  386.         
  387.     for(i=0; i<theNet->OutLayer; i++)
  388.     {    totparms += (theNet->Units[i+1])*(theNet->Units[i]);
  389.         
  390.         AllotDTypeVector(&Alpha[i], theNet->Units[i]);
  391.                     /* values for Unit[i] in each layer must be set prior
  392.                         to execution of this step.  Also note that
  393.                         the Alpha[0] vector is not used, instead a pseudo vector is
  394.                         created from a row of the XData matrix.  Thus, the allocation
  395.                         for i=0 wastes (Alpha[0].row)*sizeof(DataType) bytes.
  396.                     */
  397.         AllotDTypeVector(&dSquash[i], Alpha[i].rows);
  398.         if(Alpha[i].rows > mxprms) mxprms =  Alpha[i].rows;
  399.     }
  400.     AllotDTypeVector(&(Alpha[theNet->OutLayer]), 1);
  401.             /* this allots space for output layer, a scaler but aesthetically pleasing */
  402.             
  403.     if (theNet->method == qGaussNew)
  404.     {    AllotDTypeMatrix(&Jac_T, totparms, XData.rows+totparms);    /* actually is the transpose of Jac */
  405.         AllotDTypeVector(&Resid,XData.rows+totparms);
  406.     }
  407.     else
  408.     {    AllotDTypeMatrix(&Jac_T,totparms,XData.rows);
  409.         AllotDTypeVector(&Resid,XData.rows);
  410.     }
  411.     AllotDTypeVector(&Pi,totparms);
  412.     AllotDTypeVector(&Diag,totparms);
  413.     AllotDTypeVector(&gradSS,totparms);
  414.     AllotDTypeMatrix(&Phi,1,mxprms);    /* Phi and T2 are always a row vector since output layer is a scaler */
  415.     AllotDTypeMatrix(&T2,1,mxprms);        /* spcify as a matrix instead of vector */
  416. }
  417.  
  418. LockGradientWorkSpace()
  419. {
  420.     HLock(Resid.cells);
  421.     HLock(Jac_T.cells);
  422.     HLock(Pi.cells);
  423.     HLock(Diag.cells);
  424.     HLock(gradSS.cells);
  425.     HLockNet();
  426. }
  427.  
  428. UnlockGradientWorkSpace()
  429. {
  430.     HUnlock(Resid.cells);
  431.     HUnlock(Jac_T.cells);
  432.     HUnlock(Pi.cells);
  433.     HUnlock(Diag.cells);
  434.     HUnlock(gradSS.cells);
  435.     HUnlockNet();
  436. }